misclass_a1xz_a2xz <- function (form,y1=y1,y2=y2,   print.summary = TRUE, 
                             data = NULL) # form= y~x1+x2+x3+...+xp+z1+z2
{
  #probit <- glm(form, family = binomial(link = "probit"), data = data)
  xzmat <- model.matrix(form, data = data) #Design matrix (1,x1+x2+x3+...xp,z1,z2)
  #bmat <- probit$coef
  ysum <- model.frame(form, data = data)[, 1]
  #probitX <- glm(ysum~ xzmat[,2:(nxz-2)], family = binomial(link = "probit"), data = data)  # fit probit  considering just  x1+x2+x3+...xp
  nxz = ncol(xzmat) 
  xmat<- model.matrix(ysum~ xzmat[,2:(nxz-2)]) # Design matrix (1,x1+x2+x3+...xp)
  nx<-ncol(xmat)  
  
  # y <- model.frame(form, data = data)[, 1]
  X0<-xzmat[,1:nx]
  X1<-xzmat[,1:(nx+1)]
  X2<-cbind(xzmat[,1:nx],xzmat[,nx+2])
  
  probit.lf <- function(beta) {
    
    F0 <- pnorm(X0%*%beta[1:nx]) 
    F1 <- pnorm(X1%*%beta[(nx+1):(2*nx+1)]) 
    F2 <- pnorm(X2%*%beta[(2*nx+2):(3*nx+2)]) 
    
    
    y00 <- (1 - y1)*(1 - y2)
    y01 <- (1 - y1)*y2
    y10 <- y1*(1 - y2) 
    y11<- y1*y2
    
    
    prob00 <- 1 - F0 + F1*F2*F0 
    logprob00 <- log(prob00)
    
    prob01 <- F1*(1-F2)*F0
    logprob01 <- log(prob01)
    
    prob10 <-  (1- F1)*F2*F0
    logprob10 <- log(prob10)
    
    prob11 <-  (1- F1)*(1-F2)*F0
    logprob11 <- log(prob11)
    
    y00t <- t(y00)
    y01t <- t(y01)
    y10t <- t(y10)
    y11t <- t(y11)
    
    logl <- -sum(y00t%*%logprob00 + y01t%*%logprob01 + y10t%*%logprob10 + y11t%*%logprob11)
    return(logl)
  }
  
  
  K <- as.numeric(ncol(X0) + ncol(X1) + ncol(X2) )
  startv = rep(0,K)
  
  probitmodel.Fa1.Fa2 <- optim(startv, probit.lf, gr=NULL, method="BFGS",control=list(maxit=500), hessian=TRUE) # the function probit.lf is defined above
  coeffs.probit.Fa1.Fa2  <- probitmodel.Fa1.Fa2$par
  covmat.probit.Fa1.Fa2 <- solve(probitmodel.Fa1.Fa2$hessian)
  stderr.probit.Fa1.Fa2  <- sqrt(diag(covmat.probit.Fa1.Fa2))
  
  nameb=rep("beta",nx)
  namee1=rep("eta1",(nx+1))
  namee2=rep("eta2",(nx+1))
  
  namebeta=outer(nameb, seq(0,(nx-1)), function(x,y) paste(x,y,sep=""))[1,]
  nameeta1=outer(namee1, seq(0,(nx)), function(x,y) paste(x,y,sep=""))[1,]
  nameeta2=outer(namee2, seq(0,(nx)), function(x,y) paste(x,y,sep=""))[1,]
  
  namepar<-c(namebeta,nameeta1,nameeta2) 
  names(coeffs.probit.Fa1.Fa2 )<-namepar
  
  return(rbind(coeffs.probit.Fa1.Fa2,stderr.probit.Fa1.Fa2 ))
  
}

